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We compute the Lyapunov exponents and the Kolmogorov-Sinai (KS) entropy for a self-bound A^- 
body system that is realized as a convex billiard. This system exhibits truly high-dimensional chaos, 
and 2N — 4 Lyapunov exponents are found to be positive. The KS entropy increases linearly with 
the numbers of particles. We examine the chaos generating defocusing mechanism and investigate 
how high-dimensional chaos develops in this system with no dispersing elements. 
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I. INTRODUCTION 



^ ' Billiards are simple yet nontrivial examples of systems that display classically chaotic motion. Of special importance 
C/^ ' are the Sinai billiard |jl|] and the Bunimovich stadium j|] since they are known to be completely chaotic. Interestingly, 
these two systems exhibit two different mechanisms that generate chaos. While dispersion is the chaos-generating 
mechanism in the Sinai billiard it is defocusing that leads to chaotic dynamics in the Bunimovich stadium. Dispersing 
T-H . yields a permanent divergence of neighbored trajectories. Defocusing may occur upon reflection at a focusing boundary 
J> ' element. Provided the free path is sufficiently long nearby trajectories start to diverge after passing through the 
, focusing point, and on average the divergence might exceed the convergence thus leading to exponential instability. 
Dispersing billiards are well known also in higher dimensions. Popular examples are the three-dimensional Sinai billiard 
and the hard sphere gas. However, it was not until recently that completely chaotic billiards were constructed in more 
than two spatial dimensions that rely entirely on the defocusing mechanism These billiards use spherical caps 

as the focusing elements of the boundary. A trajectory diverges mainly in a two-dimensional plane that is defined 
. by the points of consecutive reflections with the spherical cap, and focusing may be very weak in the transversal 
■ directions. This makes it more difficult to create truly high-dimensional chaos in focusing billiards than in dispersing 
ones. Sufficient conditions for the construction of high-dimensional focusing billiards were given in ref. but it was 
1^ found that these are not necessary ones j|] . 

Besides of their intrinsic interest billiards are important model systems in the field of quantum chaos and 
statistical and fluid mechanics B. Questions related to chaos, ergodicity, transport and equilibration are often studied 
fH ' in billiard models, see e.g. refs. j^jSj. While two-dymensional chaos is fairly well understood by now, much less is known 
^ , in high-dimensional systems. Recently, a high-dimensional billiard model has been proposed in the context of nuclear 
' physics j|] and quantum chaos Within this model a self-bound A^-body system is realized as a convex billiard. 
Numerical computations yielded a positive largest Lyapunov exponent and showed that this system is predominantly 
chaotic. In this article we want to extend previous calculations and compute the full Lyapunov spectrum and the 
5^ \ KS entropy for this chaotic A^-body system. These quantities characterize the degree of hyperbolic instability in 
dynamical systems and may be related to transport coefficients in non-equilibrium situations [|ll|| . Since the studied 
billiard is convex, defocusing is the only possible source causing this instability p^ . This makes it interesting to 
examine this mechanism in more detail and compare to the situation of defocusing billiards with spherical caps. The 
results of such an investigation are not only of theoretical interest but may also be useful for further applications. We 
have in mind general questions concerning chaos in self-bound many-body systems like nuclei or atomic clusters and 
its infiuence on equilibration, damping or transport processes. 

This paper is organized as follows. In the next section we describe the model system and the techniques used to 
compute the Lyapunov exponents. The third section contains the results of our numerical computations for various 
system sizes N . In section four we investigate the defocusing mechanism in more detail. We finally give a summary. 

II. HIGH-DIMENSIONAL BILLIARD AND LYAPUNOV EXPONENTS 

Let us consider a classical system of N particles with Hamiltonian 

N 2 
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where is a two-dimensional position vector of the i-th particle and pi is its conjugate momentum. The interaction 
is given by 

rO forr<a, 
^ ' \ oo tor r > a. 

Thus, the particles move freely and interact whenever the distance between a pair of particles reaches its maximum 
value a. Hamiltonian (|l|) defines a self-bound, interacting many-body system. Energy, total momentum and total 
angular momentum are conserved quantities. For large numbers of particles the points of interactions are close to a 
circle of diameter a and therefore define a rather thin surface. Therefore, this system is a simple classical model for 
nuclei or atomic clusters. For finite values of the binding potential the system is amenable to a mean field description 
. Hamiltonian (|l|) may also be viewed as a special case of the square well gas with infinite binding potential. 
However, to the best of our knowledge, the square well gas has not been investigated for such parameter values. In 
what follows we restrict ourselves to the case of vanishing total momentum and angular momentum. 

In the limit N ^ oo the number density diverges for the self-bound many-body system A constant density 

may be obtained once the parameter a is rescaled as a ^ aN^^^, thus turning the Hamiltonian (QJ^) into an effective 
Hamiltonian. In what follows we work with a A^-independent parameter a. Since the billiard is a scaling system one 
may easily rescale the results obtained below to adapt for different values of a. 

The time evolution of a many-body system with billiard like interactions requires an effort 0{N In N) to be compared 
with the effort O^N"^) for a generic two-body interaction |jl^. Initially one computes the TV (TV — l)/2 times at which 
pairs of particles may interact and organizes these in a partially ordered binary tree, keeping the shortest time at its 
root. Immediately after an interaction of particles labelled i and j one has to recompute 2N — 3 times corresponding 
to future interactions between particles i and j and the remaining ones. The insertion of each new time into the 
partially ordered tree requires only an effort 0{\nN). Between consecutive interactions particles move freely. Upon 
an interaction of particles labelled by i and j, respectively the momenta change accordingly to 

p =p- (3) 

Here, jf and p = pi —pj are the relative momentum vectors immediately before and after the interaction, respectively, 
and r = ri — Pj is the relative position vector with magnitude |f| = a at the interaction. Obviously, eq. (||) describes 
a reflection in the center of mass system of the two interacting particles. 

We now turn to the computation of the Lyapunov exponents. We describe the used techniques rather briefly since a 
large body of literature exists on the subject, see e.g. refs. HJl^;!!^. A Hamiltonian system with / degrees of freedom 
possesses / independent Lyapunov exponents Ai, . . . , A/ ordered such that < Ai < . . . < -^/- Since the Hamiltonian 
flow preserves phase space volume there are also / non-positive Lyapunov exponents with A_j = ~Aj. A system with 
n integrals of motion has n vanishing Lyapunov exponents Ai = . . . = A„ = 0, while a chaotic system has a positive 
largest Lyapunov exponent A/ > 0. This exponent A/ is the rate at which neighbored trajectories diverge under the 
time evolution. 

Benettin et al. JTst gave a method to compute the largest Lyapunov exponent from following the time evolution of 
a reference trajectory and a second one that is initially slightly displaced. The displacement vector has to be rescaled 
after some finite evolution in a compact phase space. To compute the full spectrum of Lyapunov exponents one has 
to follow / trajectories besides the reference trajectory This defines / independent displacement vectors, and 
finite numerical precision requires their reorthogonalization besides the rescaling during the time evolution. 

Rather than following the time evolution of finite displacement vectors one may also use infinitesimal displacements 
(tangent vectors) in the computation of the Lyapunov exponents. In tangent space the time evolution is given by a 
linear mapping. Details about the tangent map in high-dimensional billiards can be found in refs. [f7|,po[. 

In a completely chaotic system the KS entropy is given by the sum of all positive Lyapunov exponents pH , i.e. 
^KS — ^j- The KS entropy measures at which rate information about the initial state of a system is lost. 

III. RESULTS 

In what follows we consider the A^-body system at vanishing total momentum and angular momentum. We use 
units such that a = m — E/N — 1. Times are then given in units of a{mN / EY/"^ . We choose initial conditions 
at random and follow a trajectory for at least 10^ collisions. This ensures a good convergence of the numerically 
computed Lyapunov spectra. 

We have checked our results as follows: The time evolution was checked by comparing forward with backward 
propagation; the Lyapunov spectra were checked by comparing the results obtained from the tangent map with those 



2 



obtained by Benettin's method involving finite displacements; the computation of all Lyapunov exponents showed that 
\-j + \j vanishes within our numerical accuracy; we found four pairs of vanishing Lyapunov exponents corresponding 
to the conserved quantities. 

The Lyapunov spectra for systems of sizes N = 10, 30, 100, 300 particles are plotted in Fig. 0. We note that the 
A^-body system possesses 2N — 4 positive Lyapunov exponents. This shows that there are no further integrals of 
motion besides energy, momentum and angular momentum, and that truly high-dimensional chaos is developed. We 
discuss this finding in detail in the following section. The Lyapunov exponent Ai is a smooth function of its index with 
a rather small smallest positive Lyapunov exponent A5 . This behavior is similar to the case of the Lennard- Jones fluid 
p2[ or the Fermi-Pasta-Ulam model |Q but differs from the hard sphere gas where a rather large smallest positive 
Lyapunov exponent was found [Q. Note that the spectra seem to converge somehow with increasing N . Table || 
displays the largest and smallest positive Lyapunov exponents, collision rates and the KS entropies. 

It is interesting to examine the iV-dependence in more detail. Fig. |^ shows that the KS entropy /iks ^^nd the 
collision rate depend linearly on the system size N . The case of the collision rate is easily understood since 
the constant single particle energy keeps the collision rate of each particle with the surface constant, too. The KS 
entropy is roughly given by the area under the corresponding spectrum presented in Fig. 0. Since the spectra converge 
approximately with increasing N this area increases linearly with the number of particles. The A^-dependence of the 
largest Lyapunov exponent A2JV is shown in Fig. ^ and may be approximated by an logarithmicly increasing curve. In 
the case of the hard sphere gas the iV-dependence could be understood for sufficiently low densities using techniques 
borrowed from kinetic theory |24|. Unfortunately, these ideas can not directly be transferred to our system since the 
density is not a small parameter. Note however, that the largest Lyapunov exponent decreases with increasing N 
once the density is kept constant after rescaling a aN^^^. This is interesting with view on nuclear physics since 
this result differs qualitatively from simple billiard (mean-field) models. Scaling arguments for such models show that 
the largest Lyapunov exponent increases with N at constant density and single-particle energy. 

The numerical results obtained in this work indicate that the considered billiard systems exhibit truly high- 
dimensional chaos. We recall that the system is convex and does not possess any dispersing elements. Furthermore, it 
differs in construction from the high-dimensional focusing billiards with spherical caps studied in refs. Thus, a 

closer examination of the chaos generating defocusing mechanism is of interest and presented in the following section. 



IV. DEFOCUSING MECHANISM 



Let us examine the defocusing mechanism in the billiard considered in this work. We do not try to proof that 
the considered system is completely chaotic - which seems difficult at least - but rather want to understand the 
numerically observed phenomenon of chaotic motion in more detail. To this purpose and based on our numerical 
results we assume that the system is (predominantly) chaotic, and that chaos is generated by the only possible 
mechanism, namely defocusing [l^ ]. We may then clarify how high-dimensional chaos develops and thus understand 
why we observe 2N — 4 positive Lyapunov exponents. This investigation may hopefully serve also as a starting point 
and a motivation for further research. 

For simplicity let us consider the three-body system first. It is useful to study this system as a billiard in full six- 
dimensional configuration space. This is possible since the change in relative momentum (|^) caused by an interaction 
of two particles corresponds to a specular reflection in the billiard. We denote vectors in conflguration space by capital 
letters as i? = {ri,r2,'r3), where = {xi,yi) is the two-dimensional position vector of the i**^ particle. The part of 
the boundary where particles labelled i = 1, 2 interact may be parametrized as 

Cq = (cosa, sina). (4) 

The (outwards pointing) normal vector daX(^i2) and the tangent vector 9q,X(i2) span the two-dimensional planes 
where divergence due to defocusing might be generated. These planes come in a four-dimensional family due to the 
parameters r and in eq. (p. Basis vectors for these planes may be chosen as 

^i = ((l,0),(-l,0),(0,0))/V2, 

^2 = ((0,1),(0,-1),(0,0))/V2. (5) 

Similar arguments show that there are two further planes where defocusing might be generated corresponding to 
interactions between particles (1,3) and (2,3), respectively. These planes are spanned by the basis vectors 
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^3 = ((l,0),(0,0),(-l,0))/\/2, 

J?4 = ((0,l),(0,0),(0,-1))/V2 (6) 



and 

^5 = ((0,0),(1,0),(-1,0))/V2, 

i6 = ((0,0),(0,l),(0,-l))/V2, (7) 

respectively. Four of the six basis vectors Ei are linearly independent. The vectors X = ((1, 0), (1, 0), (1, 0)) and 
Y = ((0, 1), (0, 1), (0, 1)) correspond to displacements of the center of mass and are orthogonal to the vectors Ei. This 
is expected since the center of mass motion moves freely. It is important to note that the boundary is neutral (i.e. 
neither focusing nor dispersing) in the transverse directions. 

It is straight forward to generalize these considerations to N bodies. In the case of the iV-body billiard there are 
TV (TV — l)/2 families of two-dimensional planes where defocusing might possibly occur. These families are related by 
those permutations that involve two out of N particles, i.e. transpositions. 2{N — 1) out of the N{N — 1) basis vectors 
Ei are linearly independent. The two vectors corresponding to the displacement of the center of mass are orthogonal 
to the vectors Ei. 

It would be interesting to relate the number of positive Lyapunov exponents and the number of linearly independent 
basis vectors Ei. Clearly, the former cannot exceed the latter. Assume that defocusing causes divergence in the 
directions of all linearly independent Ei. Then there would be exactly 2(7V — 1) positive Lyapunov exponents. 
However, the conservation of energy and angular momentum puts two additional constraints, and 2N — 4 is the 
number of positive Lyapunov exponents. This reasoning is consistent with the numerical results presented in the 
previous section. 

The following picture thus arises. The boundary of the billiard considered in this work consists of several equivalent 
elements each of which cause a reflected trajectory to diverge only in a two-dimensional plane. The orientation 
of this plane is determined by the reflecting boundary element. In transverse directions the reflection is neutral, 
i.e. neither focusing nor dispersing. A trajectory that gets reflected from sufficiently many different boundary 
elements may exhibit divergence in all directions. It is interesting to note that this mechanism differs from the one 
investigated by Bunimovich et al. |3|^. The neutral behavior in the transverse directions has the advantage that 
it avoids the problems caused by the weak convergence occurring in the transversal directions upon reflections from 
higher-dimensional spherical caps. It has the disadvantage that several focusing elements are needed to produce 
high-dimensional chaos while a single spherical cap may be sufficient. 

V. CONCLUSIONS 

We have computed the Lyapunov spectrum and the KS entropy for an interacting A^-body system in two spatial 
dimensions which is realized as a convex billiard in 2A'^-dimensional configuration space. In presence of four conserved 
quantities we find the maximal number of 2A'^ — 4 positive Lyapunov exponents. Thus, the system exhibits high- 
dimensional chaos. At fixed single particle energy the largest Lyapunov exponent grows with InA^ while the KS 
entropy grows and the collision rate increase linearly with N . In an attempt to understand the chaotic nature of 
the billiard we have identified several symmetry related two-dimensional planes where defocusing might be generated. 
Their number and orientation in configuration space is such that a long trajectory may exhibit divergence in 2^^ — 4 
directions of phase space. This mechanism of focusing differs from the one proposed recently by Bunimovich and 
Rehacek. 

Let us finally comment on chaos in realistic many-body systems. Though the considered model is a crude approx- 
imation of realistic self-bound many-body systems like nuclei or clusters it incorporates the important ingredient of 
an attractive two-body interaction that acts mainly at the surface of the system. This is the basic picture we have 
for nuclei and clusters, where the complicated two-body force creates a rather fiat mean-field potential, and particles 
experience mainly a surface interaction. However, unlike in the model system, an interaction at the nuclear surface 
involves more than just two nucleons, and the two-dimensional planes where defocusing is generated in the model 
system are replaced by some higher-dimensional ones. This might also introduce the problem caused by weak focusing 
in transverse directions. It is fair to assume that truly high-dimensional chaos may develop upon several collisions 
with the surface. Though the detailed analysis seems much more complicated than in the studied model system, the 
basic picture developed in this work should be applicable to some extend also in the case of more realistic two-body 
interactions. 
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N 


As 


A2JV 




hKS 


3 


0.70 


0.96 


2.60 


1.66 


10 


0.22 


1.14 


10.6 


11.5 


30 


0.083 


1.33 


33.1 


40.2 


100 


0.025 


1.51 


112 


141 


300 


0.008 


1.65 


338 


428 



TABLE L Smallest and largest Lyapunov exponents, collision rates and KS entropy for diflerent system sizes at fixed single 
particle energy E/N = 1. All quantities are given in units of {E /Nma^Y^'^ . 
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FIG. 1. Lyapunov spectra for N = IQ (squares), N = 30 (diamonds), A'^ = 100 (circles) and A'^ = 300 (triangles) in units of 
{E/Nma'^y^'\ The last three spectra are shifted by 0.1, 0.2 and 0.3{E/Nma'^y^^ , respectively. 



10^ e 



10^ 



00 



10^ 




10° 
10 



■ ■■! 



10' 



10' 



10" 



FIG. 2. KS entropy hsK (squares) and collision rate 1/r (circles) in units of {E/Nma ) ' as a function of system size A''. 
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